This RMarkdown file includes all code necessary to reproduce the analyses in Dablander (2025). Climate hazard experience linked to increased climate risk perception worldwide.
You can find an online data exploration tool at the World Risk Data Poll website here. We load relevant packages and the data.
library(brms)
library(dplyr)
library(haven)
library(tidyr)
library(readxl)
library(ggplot2)
library(stringr)
library(corrplot)
library(gridExtra)
library(tidyverse)
library(kableExtra)
library(countrycode)
library(RColorBrewer)
library(marginaleffects)
options(brms.backend = 'cmdstanr')
source('helpers.R')
map_hazards <- c(
'Flooding' = 'Flood\nHeavy rain\nRainstorm',
'Hurricane' = 'Hurricane\nCyclone\nTyphoon',
'Tornado' = 'Tornado',
'Thunder' = 'Thunder\nor\nlightning storm',
'Tsunami' = 'Tsunami',
'Mudslide' = 'Mudslide\nLandslide',
'Earthquake' = 'Earthquake',
'Wildfire' = 'Wildfire',
'Volcano' = 'Volcano\neruption',
'Snowstorm' = 'Blizzard\nor\nsnowstorm',
'Drought' = 'Drought',
'Heatwave' = 'Heatwave',
'Sandstorm' = 'Sandstorm'
)
map_hazards_country <- c(
'Flooding' = 'Flood | Heavy rain | Rainstorm',
'Hurricane' = 'Hurricane | Cyclone | Typhoon',
'Tornado' = 'Tornado',
'Thunder' = 'Thunder or lightning storm',
'Tsunami' = 'Tsunami',
'Mudslide' = 'Mudslide | Landslide',
'Earthquake' = 'Earthquake',
'Wildfire' = 'Wildfire',
'Volcano' = 'Volcano eruption',
'Snowstorm' = 'Blizzard or snowstorm',
'Drought' = 'Drought',
'Heatwave' = 'Heatwave',
'Sandstorm' = 'Sandstorm'
)
# We rename the relevant variables
df_full <- read.csv('data/23_wrp.csv') %>%
rename(
climate_threat = WP20719,
weather_threat = WP20723,
experienced_weather_event = WP22445,
experienced_disaster = WP23344, # Experienced a Disaster in Past Five Years in This Country
disaster_type = WP22247, # Type of Disaster Experienced in Past Five Years (Coded)
received_warning_internet = WP22248, # Received Warning About Disaster From Internet/Social Media
received_warning_government = WP22249, # Received Warning About Disaster From Local Government or Police
received_warning_media = WP22250, # Received Warning About Disaster From Radio, TV, or Newspapers
received_warning_community = WP22251, # Received Warning About Disaster From Local Community Organization
protect_future_disaster = WP22252, # Could Protect Yourself or Family in a Future Disaster
household_cover_needs = WP22228, # How Long Household Could Cover Basic Needs If Income Was Lost
household_cover_needs_weeks = WP22229, # Household Could Cover Basic Needs If Income Was Lost: Weeks
household_cover_needs_months = WP22230, # Household Could Cover Basic Needs If Income Was Lost: Months
government_cares = WP22231, # Government of (Country) Cares About You and Your Wellbeing
neighbours_care = WP22232, # Neighbors Care About You and Your Wellbeing
id = WPID_RANDOM,
country = Country,
region = GlobalRegion,
country_income = CountryIncomeLevel2023,
gender = Gender,
age = Age,
education = Education,
urbanicity = Urbanicity,
household_size = HouseholdSize,
nr_children_in_household = ChildrenInHousehold,
individual_income = INCOME_5,
resilience_individual = resilience_idv,
resilience_household = resilience_hhl,
resilience_community = resilience_com,
resilience_society = resilience_soc
) %>%
mutate(
received_warning = as.numeric(
(received_warning_internet == 1 | received_warning_government == 1 |
received_warning_media == 1 | received_warning_community == 1)
),
flooding = code_disaster(disaster_type, 1),
hurricane = code_disaster(disaster_type, 2),
wildfire = code_disaster(disaster_type, 8),
heatwave = code_disaster(disaster_type, 51),
drought = code_disaster(disaster_type, 50),
thunder = code_disaster(disaster_type, 4),
snowstorm = code_disaster(disaster_type, 10),
sandstorm = code_disaster(disaster_type, 52),
mudslide = code_disaster(disaster_type, 6),
tornado = code_disaster(disaster_type, 3),
tsunami = code_disaster(disaster_type, 5),
volcano = code_disaster(disaster_type, 9),
earthquake = code_disaster(disaster_type, 7),
experienced_disaster = recode_01(experienced_disaster),
experienced_climate_disaster = as.numeric(
flooding | hurricane | wildfire | heatwave | drought |
thunder | sandstorm | snowstorm | mudslide | tornado |
tsunami | volcano | earthquake
),
climate_threat = case_when(
climate_threat == 1 ~ 3,
climate_threat == 2 ~ 2,
climate_threat == 3 ~ 1,
TRUE ~ climate_threat
)
) %>%
select(
# Socio-demographics and similar
id, country, region, country_income, gender, age, education,
urbanicity, household_size, nr_children_in_household, individual_income,
# Threat variables and similar
climate_threat, weather_threat, experienced_weather_event, experienced_disaster,
experienced_climate_disaster, disaster_type, flooding,
hurricane, tornado, tsunami, mudslide, wildfire, drought, heatwave,
thunder, earthquake, volcano, sandstorm, snowstorm,
received_warning,
# Resilience measures and similar
household_cover_needs, household_cover_needs_weeks, household_cover_needs_months,
resilience_index, resilience_individual, resilience_household, resilience_community, resilience_society,
government_cares, neighbours_care, protect_future_disaster
)
df <- df_full
nrow(df_full)
## [1] 146910
We remove people who said that they don’t know about climate threat (coded as 98) or don’t want to answer the question (coded as 99).
table(df$climate_threat)
##
## 1 2 3 98 99
## 18804 42250 68264 17128 464
This removes about 12% of the sample.
sum(!(df$climate_threat %in% c(1, 2, 3)))
## [1] 17592
round(mean(!(df$climate_threat %in% c(1, 2, 3))), 3) * 100
## [1] 12
df <- df %>% filter(climate_threat %in% c(1, 2, 3))
We further remove people who don’t know or refuse to answer about whether they experienced a disaster in the last five years, which is 0.40%.
table(df$experienced_disaster)
##
## 0 1 98 99
## 91652 37138 485 43
sum(!(df$experienced_disaster %in% c(0, 1)))
## [1] 528
round(mean(!(df$experienced_disaster %in% c(0, 1))), 3) * 100
## [1] 0.4
df <- df %>% filter(experienced_disaster %in% c(0, 1))
We further remove people who don’t know / refuse to answer about their education (0.30%) and age (0.20%).
sum(df$education %in% c(4, 5))
## [1] 418
round(mean(df$education %in% c(4, 5)), 3) * 100
## [1] 0.3
df <- df %>% filter(!(education %in% c(4, 5)))
sum(df$age == 100)
## [1] 279
round(mean(df$age == 100), 3) * 100
## [1] 0.2
df <- df %>% filter(age != 100) # Indicates refusal
nrow(df)
## [1] 128093
The median sample size per country is \(n = 898\), with most observations from China (\(n = 2,835\)) and the least from Iceland (\(n = 458\)).
df_country <- df %>%
group_by(country) %>%
summarize(n = n()) %>%
arrange(-n)
rbind(
df_country[1, ],
df_country[142, ]
)
## # A tibble: 2 × 2
## country n
## <chr> <int>
## 1 China 2835
## 2 Iceland 458
median(df_country$n)
## [1] 898
Prepare things for the modelling further below.
df <- df %>%
mutate(
gender = factor(gender, labels = c('male', 'female')),
age = scale(age),
education = education - 1,
education = factor(education, labels = c('elementary', 'secondary', 'university')),
individual_income = factor(
individual_income,
labels = c('poorest20', 'second20', 'middle20', 'fourth20', 'richest20')
),
# Venezuela is a lower middle income country
# https://publications.iadb.org/en/venezuela-still-upper-middle-income-country-estimating-gni-capita-2015-2021?utm_source=chatgpt.com
country_income = ifelse(country_income == 9, 2, country_income),
country_income = factor(
country_income,
labels = c('lower', 'lower_middle', 'upper_middle', 'high')
)
)
Add continents.
df <- df %>%
mutate(
continent = countrycode(
country, "country.name", "continent",
custom_match = c('Kosovo' = 'Europe')
),
# Recode Americas to North or South America
continent = case_when(
country %in% c(
"United States", "Canada", "Mexico", "Cuba", "Jamaica", "Guatemala",
"El Salvador", "Honduras", "Panama", "Costa Rica", "Dominican Republic",
"Nicaragua", "Haiti", "Trinidad and Tobago", "Bahamas", "Barbados", "Belize"
) ~ "North America",
country %in% c(
"Brazil", "Argentina", "Colombia", "Chile", "Peru", "Ecuador", "Paraguay",
"Uruguay", "Venezuela", "Bolivia", "Guyana", "Suriname"
) ~ "South America",
TRUE ~ continent # Leave the rest unchanged (Europe, Africa, Asia, etc.)
)
)
df_agg <- df %>%
group_by(country) %>%
summarize(
prop_disasters = 100 * mean(experienced_disaster == 1),
prop_climate_disasters = 100 * mean(experienced_climate_disaster == 1),
prop_other_disasters = prop_disasters - prop_climate_disasters
) %>%
arrange(prop_disasters) %>%
pivot_longer(cols = -country)
df_order <- df_agg %>%
filter(name == 'prop_climate_disasters')
df_agg$country <- factor(df_agg$country, levels = df_order$country)
df_agg <- filter(df_agg, name != 'prop_disasters')
# fdablan@omics-h0.science.uva.nl:/zfs/omics/personal/fdablan
We visualize the total number of hazard experiences across all respondents and the number of countries that experienced a particular hazard.
climate_events <- c(
'flooding', 'hurricane', 'wildfire', 'heatwave',
'drought', 'thunder', 'snowstorm', 'sandstorm', 'mudslide',
'tsunami', 'volcano', 'earthquake', 'tornado'
)
df_events <- df %>%
pivot_longer(
cols = all_of(climate_events),
names_to = 'hazard',
values_to = 'value'
) %>%
group_by(country, hazard) %>%
summarise(
count = sum(value),
prop = mean(value),
.groups = "drop"
)
df_events <- df %>%
select(all_of(climate_events)) %>%
apply(2, function(x) c(sum(x), mean(x))) %>%
t() %>%
data.frame() %>%
rownames_to_column() %>%
setNames(c('event', 'count', 'proportion')) %>%
arrange(-proportion) %>%
mutate(
event = str_to_title(event),
event = map_hazards[event],
event = factor(event, levels = event)
) %>%
rename(hazard = event)
ptotal_hazards <- ggplot(df_events, aes(x = hazard, y = count, fill = 'firebrick')) +
geom_col() +
theme_minimal() +
scale_fill_manual(
values = 'firebrick'
) +
scale_y_continuous(breaks = seq(0, 16000, 2000), limits = c(0, 16000)) +
labs(x = '',y = 'Frequency') +
theme(
legend.position = 'none',
panel.grid.major = element_blank(),
panel.grid.minor = element_blank()
)
df_hazards <- df %>%
pivot_longer(
cols = all_of(climate_events),
names_to = 'hazard',
values_to = 'value'
) %>%
group_by(country, hazard) %>%
summarise(
count = sum(value),
prop = mean(value),
.groups = "drop"
) %>%
group_by(hazard) %>%
summarize(
proportion = mean(count != 0),
count = sum(count != 0)
) %>%
arrange(-count) %>%
mutate(
hazard = str_to_title(hazard),
hazard = map_hazards[hazard],
hazard = factor(hazard, levels = hazard)
)
pcountry_hazards <- ggplot(df_hazards, aes(x = hazard, y = count, fill = 'firebrick')) +
geom_col() +
theme_minimal() +
scale_y_continuous(breaks = c(seq(0, 120, 20), 142), limits = c(0, 142)) +
labs(x = '', y = 'Frequency') +
scale_fill_manual(values = 'firebrick') +
theme(
legend.position = 'none',
panel.grid.major = element_blank(),
panel.grid.minor = element_blank()
)
p2 <- ggplot(
df_hazards %>%
mutate(hazard = factor(hazard, levels = levels(df_events$hazard))),
aes(x = hazard, y = count, fill = 'firebrick')
) +
geom_col() +
theme_minimal() +
scale_y_continuous(breaks = c(seq(0, 120, 20), 142), limits = c(0, 142)) +
labs(x = '', y = 'Frequency') +
scale_fill_manual(values = 'firebrick') +
theme(
legend.position = 'none',
panel.grid.major = element_blank(),
panel.grid.minor = element_blank()
)
pdf('figures/Figure1_frequencies.pdf', width = 9, height = 9)
grid.arrange(
ncol = 1,
ptotal_hazards +
ggtitle('Total number of hazards') +
theme(plot.title = element_text(hjust = 0.50)),
p2 +
ggtitle('Number of countries with hazard') +
theme(plot.title = element_text(hjust = 0.50))
)
dev.off()
## quartz_off_screen
## 2
grid.arrange(
ncol = 1,
ptotal_hazards +
ggtitle('Total number of hazards') +
theme(plot.title = element_text(hjust = 0.50)),
p2 +
ggtitle('Number of countries with hazard') +
theme(plot.title = element_text(hjust = 0.50))
)
We visualize the percentage of people saying that climate change is a very serious threat.
tab_climate <- table(df$climate_threat)
round(tab_climate / sum(tab_climate), 3)
##
## 1 2 3
## 0.145 0.327 0.528
df_agg_threat <- df %>%
group_by(country) %>%
summarize(
prop_no = 100 * mean(climate_threat == 1),
prop_somewhat = 100 * mean(climate_threat == 2),
prop_very = 100 * mean(climate_threat == 3)
) %>%
arrange(prop_very) %>%
pivot_longer(cols = -country) %>%
mutate(
name = factor(name, levels = rev(c('prop_very', 'prop_somewhat', 'prop_no')))
)
df_order <- df_agg_threat %>%
filter(name == 'prop_very')
df_agg_threat$country <- factor(df_agg_threat$country, levels = df_order$country)
cols_world <- brewer.pal(9, 'YlOrRd')
df_agg_threat <- filter(df_agg_threat, name == 'prop_very')
breaks <- seq(0, 90, 10)
labels <- paste0('< ', breaks[seq(2, length(breaks))], '%')
df_agg_threat <- df_agg_threat %>%
mutate(
category = cut(value, breaks = breaks, labels = labels, include.lowest = TRUE),
country = as.character(country)
) %>%
remap_countries()
world_subset <- mapping_tib %>%
mutate(country_recoded = df_agg_threat$country[match(mapping_tib$country, df_agg_threat$country)]) %>%
mutate(country_recoded = ifelse(!is.na(country_recoded), country_recoded, country)) %>%
left_join(df_agg_threat, by = join_by('country'), multiple = 'first')
pworld_threat <- world_subset %>%
mutate(
text = paste('Country: ', country_recoded, '\nSample size (n): ', category, sep = '')
) %>%
ggplot(aes(x = long, y = lat, group = group, fill = factor(category), text = text)) +
geom_polygon(size = 0, alpha = 0.9) +
theme_void() +
scale_fill_manual(
values = cols_world,
na.value = '#f0f0f0', name = '',
breaks = labels,
labels = labels,
guide = guide_legend(
keyheight = unit(3, units = 'mm'),
keywidth = unit(10, units = 'mm'),
label.position = 'bottom',
title.position = 'top', nrow = 1
)) +
labs(title = 'Percentage of people considering climate change a very serious threat') +
theme(
legend.position = c(0.50, 0),
legend.title = element_text(hjust = 0.5),
plot.title = element_text(hjust = 0.5)
) +
coord_fixed(1.3)
pworld_threat
We visualize the percentage of people having a particular resilience.
df_agg_resilience <- df %>%
group_by(continent, country) %>%
summarize(resilience = mean(resilience_index, na.rm = TRUE)) %>%
arrange(-resilience) %>%
pivot_longer(cols = -c(country, continent))
df_agg_resilience$country <- factor(df_agg_resilience$country, levels = df_agg_resilience$country)
cols_world <- brewer.pal(9, 'Blues')[-seq(2)]
breaks <- seq(0.30, 0.80, 0.10)
labels <- paste0('< ', breaks[seq(2, length(breaks))])
# labels[1] <- '> 30 and < 40'
df_agg_resilience <- df_agg_resilience %>%
mutate(
category = cut(value, breaks = breaks, labels = labels, include.lowest = TRUE),
country = as.character(country)
) %>%
remap_countries()
world_subset <- mapping_tib %>%
mutate(country_recoded = df_agg_resilience$country[match(mapping_tib$country, df_agg_resilience$country)]) %>%
mutate(country_recoded = ifelse(!is.na(country_recoded), country_recoded, country)) %>%
left_join(df_agg_resilience, by = join_by('country'), multiple = 'first')
pworld_resilience <- world_subset %>%
mutate(
text = paste('Country: ', country_recoded, '\nSample size (n): ', category, sep = '')
) %>%
ggplot(aes(x = long, y = lat, group = group, fill = factor(category), text = text)) +
geom_polygon(size = 0, alpha = 0.9) +
theme_void() +
scale_fill_manual(
values = cols_world,
drop = FALSE,
na.value = '#f0f0f0', name = '',
breaks = labels,
labels = labels,
guide = guide_legend(
keyheight = unit(3, units = 'mm'),
keywidth = unit(10, units = 'mm'),
label.position = 'bottom',
title.position = 'top', nrow = 1
)) +
labs(title = 'Resilience across countries') +
theme(
legend.position = c(0.50, 0),
legend.title = element_text(hjust = 0.5),
plot.title = element_text(hjust = 0.5)
) +
coord_fixed(1.3)
pworld_resilience
Calculate resilience per continent.
df_resilience_cont <- df_agg_resilience %>%
group_by(continent) %>%
summarize(
resilience_mean = mean(value, na.rm = TRUE),
resilience_lo = quantile(value, 0.025, na.rm = TRUE),
resilience_hi = quantile(value, 0.975, na.rm = TRUE)
) %>%
arrange(-resilience_mean) %>%
mutate(
continent = factor(continent, levels = continent)
)
df_agg_resilience <- df_agg_resilience %>%
left_join(df_resilience_cont, by = 'continent') %>%
mutate(
continent = factor(continent, levels = df_resilience_cont$continent)
)
pres <- ggplot(
df_agg_resilience, aes(x = continent, y = value)
) +
geom_jitter(
aes(x = continent, y = value, color = category),
width = 0.25
) +
geom_errorbar(
aes(x = continent, ymin = resilience_lo, ymax = resilience_hi),
position = position_dodge(width = 0.25), width = 0.25, size = 1,
color = 'gray60'
) +
geom_point(
aes(x = continent, y = resilience_mean),
position = position_dodge(width = 0.25), size = 3,
color = 'gray60'
) +
# geom_col() +
# theme_minimal() +
scale_color_manual(
values = cols_world
) +
labs(x = '',y = 'Resilience') +
scale_y_continuous(breaks = seq(0.30, 0.80, 0.10), limits = c(0.30, 0.80)) +
theme_minimal() +
theme(legend.position = 'none')
pres
Combining both maps and the hazard frequencies for Figure 1.
pdf('figures/Figure1_comb.pdf', width = 18, height = 16)
grid.arrange(
ncol = 2,
pworld_threat +
theme(plot.title = element_text(hjust = 0.50, size = 18)),
ptotal_hazards +
ggtitle('Total number of hazards') +
theme(
plot.title = element_text(hjust = 0.50, size = 18),
panel.grid.major.x = element_blank(),
panel.grid.major.y = element_blank(),
panel.grid.minor.y = element_blank()
),
pworld +
theme(
legend.position = 'none',
plot.title = element_text(hjust = 0.50, size = 18)
),
p2 +
ggtitle('Number of countries with hazard') +
theme(
plot.title = element_text(hjust = 0.50, size = 18),
panel.grid.major.x = element_blank(),
panel.grid.major.y = element_blank(),
panel.grid.minor.y = element_blank()
)
)
dev.off()
## quartz_off_screen
## 2
Only resilience.
pdf('figures/Figure1a.pdf', width = 12, height = 4)
grid.arrange(
ncol = 2,
pworld_resilience +
theme(plot.title = element_text(hjust = 0.50, size = 18)),
pres +
ggtitle('Resilience across continents') +
theme(
plot.title = element_text(hjust = 0.50, size = 18)
)
)
dev.off()
## quartz_off_screen
## 2
We report the following descriptive statistics in the paper.
df_sum <- df %>%
group_by(continent) %>%
summarize(
percent_serious_threat = round(mean(climate_threat == 3), 3) * 100,
avg_climate_experience = round(mean(experienced_climate_disaster), 3) * 100,
avg_resilience = mean(resilience_index, na.rm = TRUE)
)
df_sum
## # A tibble: 6 × 4
## continent percent_serious_threat avg_climate_experience avg_resilience
## <chr> <dbl> <dbl> <dbl>
## 1 Africa 55.4 34.4 0.484
## 2 Asia 43.6 28.3 0.589
## 3 Europe 52.3 20.7 0.607
## 4 North America 64.5 22.8 0.554
## 5 Oceania 53.6 41.7 0.662
## 6 South America 73.3 24.1 0.481
We decompose the effect of a particular hazard experience into within- and between-country effects (see paper for details).
CHAINS <- 12
CORES <- 12
ITER <- 1200
WARMUP <- 500
df_between <- df %>%
group_by(country) %>%
summarize(
avg_flooding = mean(flooding),
avg_hurricane = mean(hurricane),
avg_wildfire = mean(wildfire),
avg_heatwave = mean(heatwave),
avg_drought = mean(drought),
avg_thunder = mean(thunder),
avg_snowstorm = mean(snowstorm),
avg_sandstorm = mean(sandstorm),
avg_mudslide = mean(mudslide),
avg_tsunami = mean(tsunami),
avg_volcano = mean(volcano),
avg_earthquake = mean(earthquake),
avg_tornado = mean(tornado),
avg_resilience = mean(resilience_index, na.rm = TRUE),
pct_university = mean(education == 'university')
)
df <- df %>%
left_join(df_between, by = 'country') %>%
mutate(
flooding_within = flooding - avg_flooding,
hurricane_within = hurricane - avg_hurricane,
wildfire_within = wildfire - avg_wildfire,
heatwave_within = heatwave - avg_heatwave,
drought_within = drought - avg_drought,
thunder_within = thunder - avg_thunder,
snowstorm_within = snowstorm - avg_snowstorm,
sandstorm_within = sandstorm - avg_sandstorm,
mudslide_within = mudslide - avg_mudslide,
tsunami_within = tsunami - avg_tsunami,
volcano_within = volcano - avg_volcano,
earthquake_within = earthquake - avg_earthquake,
tornado_within = tornado - avg_tornado,
resilience_within = resilience_index - avg_resilience
)
form_spec <- paste0(
'climate_threat ~ ',
'flooding_within + hurricane_within + wildfire_within + ',
'heatwave_within + drought_within + thunder_within + snowstorm_within + sandstorm_within + ',
'mudslide_within + volcano_within + earthquake_within + tornado_within + ',
'avg_flooding + avg_hurricane + avg_wildfire + ',
'avg_heatwave + avg_drought + avg_thunder + avg_snowstorm + avg_sandstorm + ',
'avg_mudslide + avg_volcano + avg_earthquake + avg_tornado + ',
'age + gender + education + individual_income + country_income + pct_university + ',
'(1 | country)'
)
df_events <- df %>% select(all_of(climate_events))
Some events are very rare, such as tsunamis. We don’t include them in the modelling.
apply(df_events, 2, sum)
## flooding hurricane wildfire heatwave drought thunder snowstorm
## 15869 4987 1681 718 3284 1065 616
## sandstorm mudslide tsunami volcano earthquake tornado
## 223 734 14 198 4992 507
round(apply(df_events, 2, mean), 3)
## flooding hurricane wildfire heatwave drought thunder snowstorm
## 0.124 0.039 0.013 0.006 0.026 0.008 0.005
## sandstorm mudslide tsunami volcano earthquake tornado
## 0.002 0.006 0.000 0.002 0.039 0.004
We run the model.
mb_spec <- run_model(
form_spec,
'models/model_multinomial.RDS',
df,
family = categorical(link = 'logit'),
chains = CHAINS, cores = CORES, iter = ITER, warmup = WARMUP
)
We compute the marginal effects and visualize them.
grid_mn <- datagrid(
model = mb_spec,
age = mean,
gender = c('male', 'female'),
education = c('elementary', 'secondary', 'university'),
country_income = unique(df$country_income),
individual_income = unique(df$individual_income)
)
df_eff_mn <- get_effects(
mb_spec, 'results/effects_main.csv',
newdata = grid_mn,
variables = c(
c(
paste0(climate_events, '_within'), paste0('avg_', climate_events),
'education', 'country_income', 'individual_income', 'pct_university'
)
), re_formula = NA
) %>%
# Only keep hazards and one contrast for education and one for income
filter(
contrast %in% c(
'dY/dX', 'university - elementary',
'richest20 - poorest20', 'high - lower', 'pct_university'
)
) %>%
mutate(
group = factor(case_when(
group == 1 ~ 'No threat at all',
group == 2 ~ 'Somewhat serious threat',
group == 3 ~ 'Very serious threat'
), levels = c('No threat at all', 'Somewhat serious threat', 'Very serious threat')),
event = ifelse(
str_starts(term, 'avg'),
str_extract(term, '(?<=_).*'),
str_extract(term, '^[^_]+')
),
climate_related = ifelse(event %in% climate_events, 1, 0),
climate_related = factor(
climate_related,
levels = c(1, 0),
labels = c('Climate-related hazard', 'Non-climate related hazard')
),
type = ifelse(str_starts(term, 'avg'), 'Between-country', 'Within-country'),
type = ifelse(term %in% c('country_income', 'pct_university'), 'Between-country', type),
type = factor(type, levels = c('Within-country', 'Between-country')),
event = str_to_title(event),
event = map_hazards[event],
event = ifelse(term == 'education', 'University\nvs\nelementary', event),
event = ifelse(term == 'pct_university', 'Percentage\nuniversity\nedcuation', event),
event = ifelse(term == 'individual_income', 'Richest 20%\nvs\npoorest 20%', event),
event = ifelse(term == 'country_income', 'High income\nvs\nlow income', event)
)
df_eff_mn_within <- df_eff_mn %>%
filter(
group == 'Very serious threat',
type == 'Within-country' # Between-country is all over the place
) %>%
arrange(climate_related, -estimate) %>%
mutate(event = factor(event, levels = event))
peff_mn <- plot_effects(df_eff_mn, fixed_levels = levels(df_eff_mn_within$event))
# Hack to use the same plotting function for the between-country effects
df_eff_mn_hack <- df_eff_mn %>%
filter(!(term %in% c('education', 'individual_income'))) %>%
mutate(
type = as.character(type),
# type = ifelse(term %in% c('education', 'individual_income'), 'Between-country', type),
estimate = ifelse(term %in% c('education', 'individual_income'), NA, estimate),
conf.low = ifelse(term %in% c('education', 'individual_income'), NA, conf.low),
conf.high = ifelse(term %in% c('education', 'individual_income'), NA, conf.high)
)
levs <- levels(df_eff_mn_within$event)
levs[c(13, 14)] <- c('Percentage\nuniversity\nedcuation', 'High income\nvs\nlow income')
peff_mn_between <- plot_effects(
df_eff_mn_hack, type = 'Between-country',
ybreaks = seq(-30, 30, 10), ylims = c(-30, 30),
fixed_levels = levs
)
grid.arrange(
ncol = 1,
peff_mn +
ggtitle('Within-country effects') +
xlab('') +
theme(
legend.position = 'none',
plot.title = element_text(hjust = 0.50),
axis.text.x = element_text(size = 8)
),
peff_mn_between +
xlab('') +
ggtitle('Between-country effects') +
theme(
legend.position = 'none',
plot.title = element_text(hjust = 0.50),
axis.text.x = element_text(size = 8)
)
)
pdf('figures/Figure2.pdf', width = 9, height = 8)
grid.arrange(
ncol = 1,
peff_mn +
ggtitle('Within-country effects') +
xlab('') +
theme(
legend.position = 'none',
plot.title = element_text(hjust = 0.50),
axis.text.x = element_text(size = 8)
),
peff_mn_between +
xlab('') +
ggtitle('Between-country effects') +
theme(
legend.position = 'none',
plot.title = element_text(hjust = 0.50),
axis.text.x = element_text(size = 8)
)
)
dev.off()
## quartz_off_screen
## 2
We run a model that includes a within-country interaction of the climate-hazard experience and an individual within-country resilience index.
df$avg_resilience_z <- as.numeric(scale(df$avg_resilience, center = TRUE, scale = TRUE))
df$resilience_within_z <- as.numeric(scale(df$resilience_within, center = TRUE, scale = TRUE))
form_spec_mod <- paste0(
'climate_threat ~ ',
'resilience_within_z * (',
'flooding_within + hurricane_within + wildfire_within + ',
'heatwave_within + drought_within + thunder_within + snowstorm_within + sandstorm_within + ',
'mudslide_within + tornado_within + volcano_within + earthquake_within) + ',
'avg_flooding + avg_hurricane + avg_wildfire + ',
'avg_heatwave + avg_drought + avg_thunder + avg_snowstorm + avg_sandstorm + ',
'avg_mudslide + avg_tornado + avg_volcano + avg_earthquake + avg_resilience_z + ',
'age + gender + education + individual_income + country_income + pct_university + ',
'(1 | country)'
)
mb_spec_mod <- run_model(
form_spec_mod,
'models/model_multinomial_moderation.RDS',
df,
family = categorical(link = 'logit'),
chains = CHAINS, cores = CORES, iter = ITER, warmup = WARMUP
)
Let’s visualize the marginal effects for different levels of resilience.
grid_mod <- datagrid(
model = mb_spec_mod,
resilience_within_z = c(-1, 0, 1),
age = mean,
gender = c('male', 'female'),
education = c('elementary', 'secondary', 'university'),
country_income = unique(df$country_income),
individual_income = unique(df$individual_income)
)
df_eff_mod <- get_effects(
mb_spec_mod, 'results/effects_moderation.csv',
newdata = grid_mod1,
variables = c(paste0(climate_events, '_within'), paste0('avg_', climate_events)),
re_formula = NA, by = 'resilience_within_z'
) %>%
# Only keep hazards and one contrast for education and one for income
mutate(
group = factor(case_when(
group == 1 ~ 'No threat at all',
group == 2 ~ 'Somewhat serious threat',
group == 3 ~ 'Very serious threat'
), levels = c('No threat at all', 'Somewhat serious threat', 'Very serious threat')),
event = ifelse(
str_starts(term, 'avg'),
str_extract(term, '(?<=_).*'),
str_extract(term, "^[^_]+")
),
type = ifelse(str_starts(term, 'avg'), 'Between-country', 'Within-country'),
type = factor(type, levels = c('Within-country', 'Between-country')),
event = str_to_title(event),
event = map_hazards[event],
resilience = ifelse(
resilience_within_z == 0, 'Average resilience',
ifelse(resilience_within_z == -1, '- 1SD resilience', '+ 1SD resilience')
),
resilience = factor(resilience, levels = c('- 1SD resilience', 'Average resilience', '+ 1SD resilience'))
)
cols <- rev(c("#D55E00", "#0072B2", "#009E73"))
pint <- ggplot(
df_eff_mod %>%
filter(
type == 'Within-country',
group == 'Very serious threat'
) %>%
mutate(
event = factor(event, levels = levels(df_eff_mn_within$event))
),
aes(x = event, y = estimate * 100, color = factor(resilience))
) +
geom_hline(yintercept = 0, linetype = "dashed") +
geom_point(position = position_dodge(width = 0.50), size = 3) +
geom_errorbar(aes(ymin = conf.low * 100, ymax = conf.high * 100),
position = position_dodge(width = 0.50), width = 0.40, size = 1) +
labs(
title = "",
x = "",
y = "Percentage difference (very serious threat)"
) +
scale_color_manual(values = cols) +
scale_y_continuous(breaks = seq(-20, 35, 5), limits = c(-20, 35)) +
theme_minimal() +
theme(
legend.position = 'top',
legend.title = element_blank(),
axis.text.x = element_text(angle = 0)
)
pint
ggsave('figures/Figure3.pdf', pint, width = 9, height = 6)
We combine the categories no threat at all and somewhat of a threat and run a logistic regression against the category very serious threat, including random slopes for countries. We do this for computational efficiency. We visualize the extent of the variation across countries and hazards.
df_country <- df %>%
mutate(
climate_threat_bin = ifelse(climate_threat == 3, 1, 0)
)
form_spec_ml <- paste0(
'climate_threat_bin ~ ',
'flooding_within + hurricane_within + wildfire_within + ',
'heatwave_within + drought_within + thunder_within + snowstorm_within + sandstorm_within + ',
'mudslide_within + volcano_within + earthquake_within + tornado_within + ',
'avg_flooding + avg_hurricane + avg_wildfire + ',
'avg_heatwave + avg_drought + avg_thunder + avg_snowstorm + avg_sandstorm + ',
'avg_mudslide + avg_volcano + avg_earthquake + avg_tornado + ',
'age + gender + education + individual_income + country_income + pct_university + ',
'(1 + flooding_within + hurricane_within + wildfire_within + ',
'heatwave_within + drought_within + thunder_within + snowstorm_within + sandstorm_within + ',
'mudslide_within + volcano_within + earthquake_within + tornado_within + ',
'education + individual_income || country)'
)
mb_country <- run_model(
form_spec_ml,
'models/model_country_slopes.RDS',
df_country,
family = bernoulli(link = 'logit'), type = 'logistic',
chains = CHAINS, cores = CORES, iter = ITER, warmup = WARMUP
)
variables <- c(paste0(climate_events, '_within'), 'individual_income', 'education')
variables <- variables[variables != 'tsunami_within']
# calculate_effects(df_country, mb_country, variables)
df_effects <- read.csv('models/effects/effects_country_percentages.csv')
df_prev <- df %>%
group_by(country) %>%
summarize(
across(
all_of(climate_events), ~ any(. == 1, na.rm = TRUE),
.names = 'has_{.col}')
) %>%
pivot_longer(
cols = starts_with('has_'),
names_to = 'hazard',
names_prefix = 'has_',
values_to = 'occurs'
)
df_effects <- df_effects %>%
mutate(
hazard = term
) %>%
left_join(
df_prev, by = c('country', 'hazard')
) %>%
mutate(
occurs = ifelse(is.na(occurs), TRUE, occurs),
climate_related = ifelse(hazard %in% c('individual', 'education'), FALSE, TRUE)
)
df_effects_group <- df_effects %>%
group_by(hazard) %>%
summarize(
estimate_mean = mean(estimate[occurs == TRUE]),
lo = quantile(estimate[occurs == TRUE], 0.025),
hi = quantile(estimate[occurs == TRUE], 0.975),
)
df_effects <- df_effects %>%
left_join(df_effects_group, by = 'hazard')
df_effects_occurs <- filter(df_effects, occurs)
df_effects_order <- df_effects_group %>%
mutate(size_ci = hi - lo) %>%
arrange(size_ci)
hazard_order <- c(
df_effects_order$hazard[!(df_effects_order$hazard %in% c('individual', 'education'))],
'education', 'individual'
)
df_effects_occurs$hazard <- factor(
df_effects_occurs$hazard,
levels = hazard_order
)
map_all <- map_hazards
map_all['Individual'] <- 'Richest 20%\nvs\npoorest 20%'
map_all['Education'] <- 'University\nvs\nelementary'
df_effects_occurs <- df_effects_occurs %>%
mutate(
hazard = factor(
hazard,
labels = map_all[str_to_title(levels(hazard))]
# labels = levs
)
)
set.seed(1)
pcountry <- ggplot(
data = df_effects_occurs,
aes(x = hazard, y = estimate * 100, col = factor(climate_related))
) +
geom_jitter(
aes(x = hazard, y = estimate * 100),
# col = 'gray90'
col = 'black', alpha = 0.25
) +
geom_point(
aes(x = hazard, y = estimate_mean * 100),
position = position_dodge(width = 0.25), size = 3
# col = 'firebrick'
) +
geom_errorbar(
aes(x = hazard, ymin = lo * 100, ymax = hi * 100),
position = position_dodge(width = 0.25), width = 0.25, size = 1,
# col = 'firebrick'
) +
labs(
title = "",
x = "",
y = "Percentage difference (very serious threat)",
fill = "Disaster Type"
) +
scale_y_continuous(breaks = seq(-10, 25, 5), limits = c(-10, 25)) +
scale_color_manual(values = c('gray40', 'firebrick')) +
theme_minimal() +
theme(
legend.position = 'none',
strip.text.x = element_text(size = 12),
axis.text.y = element_text(size = 8),
legend.title = element_blank()
# panel.grid.minor.x = element_blank(),
# panel.grid.major.x = element_blank()
)
pcountry
ggsave('figures/Figure4.pdf', pcountry, width = 11, height = 5)
The data includes a variable on whether the respondents have experienced an extreme weather event in the last two years. We analyse the effect of this on climate risk perceptions here.
df_ext <- df %>%
filter(experienced_weather_event %in% c(1, 2, 3, 4)) %>%
mutate(
ew_personally = ifelse(experienced_weather_event == 1, 1, 0),
ew_somebody = ifelse(experienced_weather_event == 2, 1, 0),
ew_both = ifelse(experienced_weather_event == 3, 1, 0),
ew_none = ifelse(experienced_weather_event == 4, 1, 0)
)
df_extc <- df_ext %>%
group_by(country) %>%
summarize(
avg_ew_none = mean(ew_none),
avg_ew_personally = mean(ew_personally),
avg_ew_somebody = mean(ew_somebody),
avg_ew_both = mean(ew_both)
)
df_ext <- df_ext %>%
left_join(df_extc, by = 'country') %>%
mutate(
ew_none_within = ew_none - avg_ew_none,
ew_personally_within = ew_personally - avg_ew_personally,
ew_somebody_within = ew_somebody - avg_ew_somebody,
ew_both_within = ew_both - avg_ew_both
)
table(df_ext$experienced_weather_event)
##
## 1 2 3 4
## 11213 22502 5674 88117
df_ext <- df %>%
# 1: Yes, personally
# 2: Yes, somebody I know
# 3: Both
# 4: No
filter(experienced_weather_event %in% c(1, 3, 4)) %>%
mutate(
experienced_extreme = ifelse(experienced_weather_event %in% c(1, 3), 1, 0)
)
df_extc <- df_ext %>%
group_by(country) %>%
summarize(avg_experienced_extreme = mean(experienced_extreme))
df_ext <- df_ext %>%
left_join(df_extc, by = 'country') %>%
mutate(experienced_extreme_within = experienced_extreme - avg_experienced_extreme)
form_extreme <- paste0(
'climate_threat ~ ',
'experienced_extreme_within + avg_experienced_extreme + ',
'age + gender + education + individual_income + country_income + pct_university + ',
'(1 | country)'
)
mb_ext <- run_model(
form_extreme,
'models/model_multinomial_ew.RDS',
df_ext,
family = categorical(link = 'logit'),
chains = CHAINS, cores = CORES, iter = ITER, warmup = WARMUP
)
grid_ew <- datagrid(
model = mb_ext,
age = mean,
gender = c('male', 'female'),
education = c('elementary', 'secondary', 'university'),
country_income = unique(df$country_income),
individual_income = unique(df$individual_income)
)
df_eff_ew <- get_effects(
mb_ext, 'results/effects_ew.csv',
newdata = grid_ew,
variables = c(
c(
'experienced_extreme_within', 'avg_experienced_extreme',
'education', 'country_income', 'individual_income', 'pct_university'
)
), re_formula = NA
)
# 11.7% difference in people considering climate change a very serious threat
df_eff_ew %>%
filter(group == 3)
## term group contrast estimate
## 1 avg_experienced_extreme 3 dY/dX 0.2706476365
## 2 country_income 3 high - lower 0.0486850481
## 3 country_income 3 lower_middle - lower 0.0146105839
## 4 country_income 3 upper_middle - lower 0.0510359081
## 5 education 3 secondary - elementary 0.0488498544
## 6 education 3 university - elementary 0.1079579607
## 7 experienced_extreme_within 3 dY/dX 0.1172912595
## 8 individual_income 3 fourth20 - poorest20 0.0009423518
## 9 individual_income 3 middle20 - poorest20 0.0002188238
## 10 individual_income 3 richest20 - poorest20 0.0100864561
## 11 individual_income 3 second20 - poorest20 0.0032031282
## 12 pct_university 3 dY/dX -0.2223871970
## conf.low conf.high
## 1 -0.1022411677 0.61985399
## 2 -0.1271273092 0.21847175
## 3 -0.1338698520 0.14968071
## 4 -0.0980840840 0.19532870
## 5 0.0401301783 0.05799999
## 6 0.0963859153 0.11956600
## 7 0.1075422651 0.12696092
## 8 -0.0095591995 0.01150181
## 9 -0.0104248461 0.01091673
## 10 -0.0002260002 0.02045212
## 11 -0.0075901392 0.01437635
## 12 -0.5389898443 0.11143318
Running a multilevel logistic regression to assess country heterogeneity.
form_extreme_bin <- paste0(
'climate_threat_bin ~ ',
'experienced_extreme_within + avg_experienced_extreme + ',
'age + gender + education + individual_income + country_income + pct_university + ',
'(1 + experienced_extreme_within + education + individual_income || country)'
)
df_ext$climate_threat_bin <- ifelse(df_ext$climate_threat == 3, 1, 0)
mb_ext <- run_model(
form_extreme_bin,
'models/model_country_slopes_ew.RDS',
df_ext,
family = bernoulli(link = 'logit'), type = 'logistic',
chains = CHAINS, cores = CORES, iter = ITER, warmup = WARMUP
)
Let’s estimate the country variation in extreme weather events.
# calculate_effects(df_ext, mb_ext, 'experienced_extreme_within')
df_effects_ew <- read.csv('models/effects/effects_country_percentages_ew.csv')
quantile(df_effects_ew$estimate, c(0.025, 0.975))
## 2.5% 97.5%
## 0.03966262 0.17982184
df_country_disasters <- df %>%
group_by(country) %>%
summarize(
flooding_prop = mean(flooding),
hurricane_prop = mean(hurricane),
wildfire_prop = mean(wildfire),
heatwave_prop = mean(heatwave),
drought_prop = mean(drought),
thunder_prop = mean(thunder),
snowstorm_prop = mean(snowstorm),
sandstorm_prop = mean(sandstorm),
mudslide_prop = mean(mudslide),
tsunami_prop = mean(tsunami),
tornado_prop = mean(tornado),
volcano_prop = mean(volcano),
earthquake_prop = mean(earthquake)
) %>%
pivot_longer(cols = -country)
map1 <- make_disaster_map(df_country_disasters, disaster = 'heatwave', title = 'Heatwave')
map2 <- make_disaster_map(df_country_disasters, disaster = 'mudslide', title = 'Mudslide | Landslide')
map3 <- make_disaster_map(df_country_disasters, disaster = 'sandstorm', title = 'Sandstorm')
map4 <- make_disaster_map(df_country_disasters, disaster = 'flooding', title = 'Flood | Heavy rain | Rainstorm')
map5 <- make_disaster_map(df_country_disasters, disaster = 'drought', title = 'Drought')
map6 <- make_disaster_map(df_country_disasters, disaster = 'thunder', title = 'Thunder or lightning storm')
map7 <- make_disaster_map(df_country_disasters, disaster = 'wildfire', title = 'Wildfire')
map8 <- make_disaster_map(df_country_disasters, disaster = 'hurricane', title = 'Hurricane | Cyclone | Typhoon')
map9 <- make_disaster_map(df_country_disasters, disaster = 'tornado', title = 'Tornado')
map10 <- make_disaster_map(df_country_disasters, disaster = 'snowstorm', title = 'Blizzard or snowstorm')
map11 <- make_disaster_map(df_country_disasters, disaster = 'earthquake', title = 'Earthquake')
map12 <- make_disaster_map(df_country_disasters, disaster = 'volcano', title = 'Volcano eruption')
pdf('figures/FigureS1_disasterall.pdf', width = 8, height = 10)
grid.arrange(
ncol = 3,
map1, map2, map3, map4, map5, map6, map7, map8, map9, map10, map11, map12
)
dev.off()
## quartz_off_screen
## 2
We visualize the within- and between-country effects across all levels of climate risk perception.
cols <- rev(c("#D55E00", "#0072B2", "#009E73"))
pallw <- ggplot(
df_eff_mn %>%
filter(type == 'Within-country') %>%
mutate(event = factor(event, levels = levels(df_eff_mn_within$event))),
aes(x = event, y = estimate * 100, color = factor(group))
) +
geom_hline(yintercept = 0, linetype = "dashed") +
geom_point(position = position_dodge(width = 0.50), size = 3) +
geom_errorbar(aes(ymin = conf.low * 100, ymax = conf.high * 100),
position = position_dodge(width = 0.50), width = 0.40, size = 1) +
labs(
title = "",
x = "",
y = "Percentage difference (very serious threat)"
) +
scale_color_manual(values = cols) +
scale_y_continuous(breaks = seq(-20, 20, 5), limits = c(-20, 20)) +
theme_minimal() +
theme(
legend.position = 'top',
legend.title = element_blank(),
axis.text.x = element_text(angle = 0)
)
levs <- levels(df_eff_mn_within$event)
levs[c(13, 14)] <- c('Percentage\nuniversity\nedcuation', 'High income\nvs\nlow income')
pallb <- ggplot(
df_eff_mn_hack %>%
filter(type == 'Between-country') %>%
mutate(event = factor(event, levels = levs)),
aes(x = event, y = estimate * 10, color = factor(group))
) +
geom_hline(yintercept = 0, linetype = "dashed") +
geom_point(position = position_dodge(width = 0.50), size = 3) +
geom_errorbar(aes(ymin = conf.low * 10, ymax = conf.high * 10),
position = position_dodge(width = 0.50), width = 0.40, size = 1) +
labs(
title = "",
x = "",
y = "Percentage difference (very serious threat)"
) +
scale_color_manual(values = cols) +
scale_y_continuous(breaks = seq(-40, 40, 10), limits = c(-40, 40)) +
theme_minimal() +
theme(
legend.position = 'top',
legend.title = element_blank(),
axis.text.x = element_text(angle = 0)
)
grid.arrange(
ncol = 1,
pallw +
ggtitle('Within-country effects') +
xlab('') +
theme(
legend.position = 'none',
plot.title = element_text(hjust = 0.50),
axis.text.x = element_text(size = 8)
),
pallb +
xlab('') +
ggtitle('Between-country effects') +
theme(
legend.position = 'none',
plot.title = element_text(hjust = 0.50),
axis.text.x = element_text(size = 8)
)
)
pdf('figures/FigureS2.pdf', width = 9, height = 10)
grid.arrange(
ncol = 1,
pallw +
ggtitle('Within-country effects') +
xlab('') +
theme(
plot.title = element_text(hjust = 0.50),
axis.text.x = element_text(size = 8)
),
pallb +
xlab('') +
ggtitle('Between-country effects') +
theme(
# legend.position = 'none',
plot.title = element_text(hjust = 0.50),
axis.text.x = element_text(size = 8)
)
)
dev.off()
## quartz_off_screen
## 2
We visualize the Pearson correlations of within- and between-country as well as socio-demographic variables.
preds <- c(
'flooding_within', 'hurricane_within', 'wildfire_within',
'heatwave_within', 'drought_within', 'thunder_within', 'snowstorm_within', 'sandstorm_within',
'mudslide_within', 'tornado_within', 'volcano_within', 'earthquake_within', 'resilience_within',
'avg_flooding', 'avg_hurricane', 'avg_wildfire', 'avg_heatwave', 'avg_drought',
'avg_thunder', 'avg_snowstorm', 'avg_sandstorm', 'avg_mudslide', 'avg_tornado',
'avg_volcano', 'avg_earthquake', 'avg_resilience', 'age', 'gender', 'education',
'individual_income', 'country_income', 'pct_university'
)
df_preds <- df %>%
select(all_of(preds)) %>%
mutate_all(as.numeric)
cormat <- cor(df_preds, method = 'pearson', use = 'pairwise.complete.obs')
diag(cormat) <- NA
cnames <- c(
'Flood | Heavy rain | Rainstorm (within)',
'Hurricane | Cyclone | Typhoon (within)',
'Wildfire (within)',
'Heatwave (within)',
'Drought (within)',
'Thunder or lightning storm (within)',
'Blizzard or snowstorm (within)',
'Sandstorm (within)',
'Mudslide | Landslide (within)',
'Tornado (within)',
'Volcano eruption (within)',
'Earthquake (within)',
'Resilience (within)',
'Flood | Heavy rain | Rainstorm (between)',
'Hurricane | Cyclone | Typhoon (between)',
'Wildfire (between)',
'Heatwave (between)',
'Drought (between)',
'Thunder or lightning storm (between)',
'Blizzard or snowstorm (between)',
'Sandstorm (between)',
'Mudslide | Landslide (between)',
'Tornado (between)',
'Volcano eruption (between)',
'Earthquake (between)',
'Resilience (between)',
'Age', 'Gender', 'Education', 'Individual income', 'Country income', '% University'
)
colnames(cormat) <- rownames(cormat) <- cnames
pdf('figures/FigureS3.pdf', width = 10, height = 9)
# Hazards trivially not correlated because one can only indicate a single hazard
corrplot(
cormat, method = 'color', type = 'upper', number.cex = 0.50,
addCoef.col = 'black',
tl.cex = 0.60, addrect = 20, tl.col = 'black',
na.label = ' '
)
dev.off()
## quartz_off_screen
## 2
corrplot(
cormat, method = 'color', type = 'upper', number.cex = 0.50,
addCoef.col = 'black',
tl.cex = 0.60, addrect = 20, tl.col = 'black',
na.label = ' '
)
peffc <- ggplot(
df_effects_occurs,
aes(x = country, y = estimate * 100, color = factor(climate_related))
) +
geom_hline(aes(yintercept = estimate_mean * 100), col = 'gray76') +
geom_point(
position = position_dodge(width = 0.25), size = 0.60
) +
geom_errorbar(
aes(ymin = conf.low * 100, ymax = conf.high * 100),
position = position_dodge(width = 0.25),
width = 0.75, size = 0.4
) +
labs(
title = "",
x = "",
y = "Percentage difference (very serious threat)",
fill = "Disaster Type"
) +
scale_color_manual(values = c('gray40', 'firebrick')) +
facet_wrap(~ hazard, ncol = 2) +
scale_y_continuous(breaks = seq(-20, 30, 10)) +
coord_cartesian(ylim = c(-20, 30)) + # <--- this is the key line
theme_minimal() +
theme(
legend.position = 'none',
strip.text.x = element_text(size = 12),
# axis.text.x = element_text(size = 2.5, angle = 90),
axis.text.x = element_blank(),
axis.text.y = element_text(size = 8),
legend.title = element_blank(),
panel.grid.minor.x = element_blank(),
panel.grid.major.x = element_blank()
)
peffc
ggsave('figures/FigureS4.pdf', peffc, width = 9, height = 14)
df_sd <- VarCorr(mb_country)$country$sd %>%
data.frame() %>%
rownames_to_column() %>%
filter(
rowname %in% c(variables, 'individual_incomerichest20', 'educationuniversity')
) %>%
mutate(
term = map_chr(str_split(rowname, '_'), 1),
term = ifelse(
!(term %in% c('individual', 'educationuniversity')),
map_hazards[str_to_title(term)], term
),
term = ifelse(term == 'individual', 'Richest 20%\nvs\npoorest 20%', term),
term = ifelse(term == 'educationuniversity', 'University\nvs\nelementary', term)
)
df_sd <- bind_rows(
df_sd %>%
filter(!(term %in% c('University\nvs\nelementary', 'Richest 20%\nvs\npoorest 20%'))) %>%
arrange(Estimate) %>%
mutate(climate_related = TRUE),
df_sd[c(13, 14), ] %>%
mutate(climate_related = FALSE)
) %>%
mutate(
term = factor(term, levels = term)
)
psd <- ggplot(
data = df_sd,
aes(x = term, y = Estimate, col = factor(climate_related))
) +
geom_point(
aes(x = term, y = Estimate),
col = 'black', alpha = 0.25
) +
geom_point(
aes(x = term, y = Estimate),
position = position_dodge(width = 0.25), size = 3
) +
geom_errorbar(
aes(x = term, ymin = Q2.5, ymax = Q97.5),
position = position_dodge(width = 0.25), width = 0.25, size = 1
) +
labs(
title = "",
x = "",
y = "Standard deviation of random slope",
fill = "Disaster Type"
) +
scale_y_continuous(breaks = seq(0, 2, 0.10), limits = c(0, 2.1)) +
coord_cartesian(ylim = c(0, 0.80)) + # <--- this is the key line
scale_color_manual(values = c('gray40', 'firebrick')) +
theme_minimal() +
theme(
legend.position = 'none',
strip.text.x = element_text(size = 12),
axis.text.y = element_text(size = 8),
legend.title = element_blank()
)
ggsave('figures/FigureS5.pdf', psd, width = 10, height = 5)
sessionInfo()
## R version 4.3.3 (2024-02-29)
## Platform: aarch64-apple-darwin20 (64-bit)
## Running under: macOS Sonoma 14.3
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.11.0
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## time zone: Europe/Stockholm
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] marginaleffects_0.25.0 RColorBrewer_1.1-3 countrycode_1.6.1
## [4] kableExtra_1.4.0 lubridate_1.9.3 forcats_1.0.0
## [7] purrr_1.0.2 readr_2.1.5 tibble_3.2.1
## [10] tidyverse_2.0.0 gridExtra_2.3 corrplot_0.92
## [13] stringr_1.5.1 ggplot2_3.5.1 readxl_1.4.3
## [16] tidyr_1.3.1 haven_2.5.4 dplyr_1.1.4
## [19] brms_2.21.0 Rcpp_1.0.12
##
## loaded via a namespace (and not attached):
## [1] tidyselect_1.2.1 viridisLite_0.4.2 farver_2.1.1
## [4] loo_2.7.0 fastmap_1.1.1 tensorA_0.36.2.1
## [7] digest_0.6.35 estimability_1.5.1 timechange_0.3.0
## [10] lifecycle_1.0.4 StanHeaders_2.32.6 processx_3.8.4
## [13] magrittr_2.0.3 posterior_1.5.0 compiler_4.3.3
## [16] rlang_1.1.5 sass_0.4.9 tools_4.3.3
## [19] utf8_1.2.4 yaml_2.3.8 data.table_1.17.0
## [22] knitr_1.45 labeling_0.4.3 bridgesampling_1.1-2
## [25] pkgbuild_1.4.4 curl_5.2.1 cmdstanr_0.9.0
## [28] xml2_1.3.6 abind_1.4-5 withr_3.0.0
## [31] grid_4.3.3 stats4_4.3.3 fansi_1.0.6
## [34] xtable_1.8-4 colorspace_2.1-0 inline_0.3.19
## [37] emmeans_1.10.2 scales_1.3.0 insight_1.1.0
## [40] cli_3.6.2 mvtnorm_1.2-4 rmarkdown_2.27
## [43] ragg_1.3.0 generics_0.1.3 RcppParallel_5.1.7
## [46] rstudioapi_0.16.0 tzdb_0.4.0 cachem_1.0.8
## [49] rstan_2.32.6 bayesplot_1.11.1 parallel_4.3.3
## [52] cellranger_1.1.0 matrixStats_1.3.0 vctrs_0.6.5
## [55] V8_4.4.2 Matrix_1.6-5 jsonlite_1.8.8
## [58] hms_1.1.3 systemfonts_1.0.6 jquerylib_0.1.4
## [61] glue_1.7.0 ps_1.7.6 codetools_0.2-19
## [64] distributional_0.4.0 stringi_1.8.3 gtable_0.3.5
## [67] QuickJSR_1.1.3 munsell_0.5.1 pillar_1.9.0
## [70] htmltools_0.5.8.1 Brobdingnag_1.2-9 R6_2.5.1
## [73] textshaping_0.3.7 evaluate_0.23 lattice_0.22-5
## [76] highr_0.10 backports_1.4.1 bslib_0.7.0
## [79] rstantools_2.4.0 svglite_2.1.3 coda_0.19-4.1
## [82] nlme_3.1-164 checkmate_2.3.1 xfun_0.43
## [85] pkgconfig_2.0.3